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A Local Representation of the Electronic Dielectric Response Function 
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We present a local representation of the electronic dielectric response function, based on a spatial 
partition of the dielectric response into contributions from each Wannier function using a generalized 
density functional perturbation theory. This procedure is fully ah initio, and therefore allows us to 
rigorously define local metrics, such as “bond polarizability”, on Wannier centers. We show that the 
locality of the response function is determined by the locality of three quantities: Wannier functions 
of the occupied manifold, the density matrix, and the Hamiltonian matrix. In systems with a gap, 
the bare dielectric response is exponentially localized, which supports the physical picture of the 
dielectric response function as a collection of interacting local response that can be captured by a 
tight-binding model. 

PACS numbers: 71.15.-m, 71.45.Gm, 71.15.Qe, 71.15.Ap 


The screened electronic dielectric response function 
(EDRF), X, is a fundamental physical quantity that cap¬ 
tures many-electron correlation effects. From a micro¬ 
scopic point of view, y relates the perturbation from an 
external potential at r' to the electronic density response 
at r. This intrinsic non-local character has precluded 
a compact local representation of y in electronic struc¬ 
ture theory; it has been traditionally represented by huge 
matrices in either real or reciprocal space [1]. This cum¬ 
bersome matrix representation of y has become a ma¬ 
jor computational bottleneck to accurately predict elec¬ 
tron correlation energy and electronic excitation spec¬ 
tra. More importantly, the physical interpretation of y 
is largely limited to its macroscopic average. A robust, 
microscopic theory that describes the local characteris¬ 
tics (e.g., shape, strength and decay rate) is needed to 
unravel the underlying physical nature of EDRFs. 

Empirical methods have been used to partition EDRFs 
to obtain polarizabilities or effective van der Waals Cq 
dispersion coefficients of atoms inside either a molecule 
or a solid [2, 3]. On the other hand, non-empirical 
methods often rely on extra approximations to parti¬ 
tion, e.g., the molecular polarizability, into so-called “dis¬ 
tributed polarizabilities” [4]. Such procedures typically 
involve partitioning the volume [5] or the basis space [6] 
of a molecule, or fitting the point-to-point polarizabil¬ 
ities computed on a grid around a molecule [7], with 
known drawbacks including large charge-flow terms that 
are hard to localize, strong basis set dependence, and 
high computational cost [8]. 

The Wannier function (WE) representation [9] is a nat¬ 
ural choice to describe chemical bonds using either “Boys 
orbitals” for molecules [10] or the maximally localized 
Wannier functions (MLWFs) for the solid-state equiva¬ 
lent [11-13]. However, the link between local EDRFs 
and WFs is obscured by the fact that EDRFs concern 
electron-hole pairs rather than electronic orbitals alone. 
Silvestrelli [14] proposed to define Cq coefficients on Wan¬ 
nier centers using empirical models, assuming that WFs 


have the 5-symmetry. Giustino and Pasquarello [15] in¬ 
troduced the local dielectric permittivity in layered sys¬ 
tems, based on local dipole moments derived from each 
Wannier function using the Berry-phase theory of the 
polarization [16]. However, this approach requires sepa¬ 
rate calculations under finite external fields, unsuitable 
to study the spatial decay rate and the dynamic response. 
Lu et al. [17] applied a simultaneous diagonalization algo¬ 
rithm to directly localize the eigenvectors of EDRFs. De¬ 
spite the observed trend in the locality [17], the chemical 
nature of localized EDRFs was not determined precisely. 

In this Letter, we propose a local representation for 
microscopic EDRFs using a generalized density func¬ 
tional perturbation theory (DEPT) [18]. While the con¬ 
ventional theory is formulated on the eigenstates of the 
Kohn-Sham (KS) Hamiltonian, we generalize the DEPT 
to any orthogonal basis set that spans the occupied state 
manifold. A convenient choice adopted in this work is 
the MLWF [11, 12], as it can provide insightful interpre¬ 
tations regarding chemical bonds. Because this method 
is fully ab initio, it ensures accuracy and transferability. 

First we generalize DEPT for the bare EDRF, y^, 
which is the building block in linear response theory. 
Under a perturbation in the self-consistent potential, 
the response density can be calculated 

through y^ as 


Ap{u};r) = j dr'x°(w;r,r') AV;(r'). (1) 

In the following, we adopt the shorthand notation: Ap = 
y^ AVg. X can be solved from y^ through Dyson’s equa¬ 
tion, y = y^ + y^ y, where K = VcS- K^c with Vc and 
Kxc being Goulomb and exchange-correlation kernels, re¬ 
spectively [19]. For periodic systems, Eq. 1 is often solved 
for individual Fourier components with wave vectors q, 
AI/^(r) = Avs{t). The response density matrix is 
given by 


§ X {^v 

vk 


( 2 ) 
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where the variation of the KS orbital, is the 

solution of the Sternheimer equation [18], 


aP'^+^±oj) I. (3) 


Here makes Eq. 3 non singular; and P^^^ 

are projectors onto the occupied and unoccupied state 
manifolds at momentum k + q, which are introduced to 
avoid the explicit reference to the unoccupied states [18]. 

It is trivial to partition Apq in Eq. 2 in the energy 
domain into contributions from individual KS orbitals at 
given {v^ k} and q, and the corresponding linear equa¬ 
tions in Eq. 3 are decoupled. Alternatively, a real space 
partition of EDREs can be achieved through a general¬ 
ized DEPT in the WE representation. Eollowing the no¬ 
tations in Ref. [11], we define WEs and their first order 
perturbations as 


WR„(r) = f dke 

(27r) Jbz 

Jbz q 

( 4 ) 

where J is the number of orbitals used to construct WEs, 
and unitary matrices minimize the spatial spreads 
of the WEs labeled by lattice vector R and band index 
n [11]. Eor simplicity, we will focus on the the static limit 
and drop the superscripts + and —; extension to the dy¬ 
namic case is straightforward. Because p and, thus, Ap 
are invariant under the unitary transformation of occu¬ 
pied states, Eq. 2 can be rewritten in the Wannier repre¬ 
sentation as 


Rn 

Applying to both sides of Eq. 3 and integrating 
over k, one obtains the generalized Sternheimer equation 
in the Wannier representation as 

E (^'Rn,R'n' -H- aa) |AW^r.„.) = PcAV |W^Rn) . 

R'n' 

(6) 

Because WEs are not eigenstates of the KS Hamiltonian, 
unlike Eq. 3, Eq. 6 yields a set of coupled equations due to 
the hopping integral terms, 6Rn,R'n' = (f^Rnl^lVER'n') 
(Rn 7^ R'n'). It follows that 


|AtyR„) = ^ \AWnn,n'n') 

R'n' 

^J 2 [e-{H + aP,) JlR^'n'^AV |W^R'n'), 

R'n' 

( 7 ) 

where I is an Nw x Nw identity matrix with Nw be¬ 
ing the number of occupied WEs. AlERn,R'n' denotes 


the variation of Wyiu caused by the perturbation at 
Hpi'n'- Combining Eqs. 5 and 7, one can expand 
in terms of two-body partial response functions, — 
ERn,R'n' Xrh.r'hm with Corresponding density response 
given by = ^|AWr„,r/„/)(Wr„|. A formal 

real space partition can be established by contracting 
Xru R'n' one-body variables, 

X°(r,r') = y]xR„(r,r'), 

Rn 

XRn(l’,l’') = \X1 [XR„,R'n'(l’,I’') + XR/„yR„(r,r')] • 

^ R'n' 

(8) 

Hermitian by construction. It contains not only 
on-site terms, but also the charge-flow (R'n' ^ Rn) into 
and out of the local orbital, whose magnitude determines 
the extent of the locality of Eq. 8 is the central 
equation of our theory, which together with the local¬ 
ity analysis below formally establishes the local repre¬ 
sentation of EDREs. The real space partition can be 
achieved similarly for y through partial response func- 
tions, XR„ = 4 X]r,„, (xRn.R'n' + XR'n',Rn), where 

CX) 

XRn,R'n' ^ ^ ^ > XRn,n"n" X )R"n",R'n' ’ 

Yl"n" m=0 

CX) 

XR'n',Rn = (x ^)R/n',n"n" X-R"n",'Rn^{^) 

Yl"n" m=0 

are solved self-consistently with Eq. 6. We emphasize 
that there is a critical distinction between Eqs. 8-9 and 
existing approximate partition methods. In our method, 
once a Wannier localization procedure is chosen, the sub¬ 
sequent real space EDRE partition is exact. 

To demonstrate the concept of the local EDRE, we 
partition static molecular polarizabilities of acetylene 
{C 2 H 2 ), ethylene {C 2 H 4 ), and ethane {C 2 HQ) into bond 
polarizabilities at Wannier centers. The major differ¬ 
ence between these molecules is the C-C bond order: 
single bond in C 2 HQ, double bond in C 2 H 4 , and triple 
bond in C 2 H 2 . We compare the standard mean screened 
polarizability {a = |tr(rxf')) with the unscreened one 
= |tr(rx^f')) to gain insight into the local field ef¬ 
fect. All the calculations were performed at the random 
phase approximation level using WannierQO [20] and a 
modified version of Quantum ESPRESSO [21]; com¬ 
putational details are given in the Supplemental Material 
(SM) [22]. 

Within C 2 H 2 , as shown in Table I, a^{7rcc) > 
a^{CH) > a^{acc)‘ This behavior is a direct outcome 
of the electronic structure, because ircc is closest to the 
lowest unoccupied molecular orbital (LUMO) in energy 
(i.e., most reactive), and acc is farthest from LUMO 
(i.e., least reactive) as shown by the projected density of 
states in Eig. S2 in SM. The same trend holds for all three 
molecules. Among different molecules, unscreened bond 
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TABLE I: Mean static molecular and bond polarizabilities 
(in Bohr^) of C 2 HQ, C 2 H 4 , and C 2 H 2 . Multipliers in bond 
polarizabilities indicate the degeneracy; unscreened bond po¬ 
larizabilities are shown in the parentheses. 



Total C-H C-C 

acc 7 ^cc 

C2H2 

C2H4 

C2HG 

24.65 2x3.11(4.36) 18.43 
29.41 4x3.91(5.73) 13.76 
31.35 6x4.54(6.60) 4.09 

0.60(1.55) 2x8.92(16.26) 
1.73(3.30) 12.02(21.28) 

4.09(6.25) 


polarizabilities increase in the order of C 2 H 2 , C' 2 i^ 4 , and 
C 2 Hq^ partially due to the reduced bonding-antibonding 
splitting with increased bond lengths {dcH' 1 . 000 , 1.018, 
1.024; dec'' 1-000, 1.102, 1.267, both normalized by the 
bond lengths of C 2 H 2 ). Another dominating factor in 
the molecular polarizability is the number of bonds, i.e., 
the degeneracy. On the other hand, a is always smaller 
than q;^, because of the screening effect. We define 
Ceff = a:^/(a as a measure of the local screening strength 
that includes both intra- and inter-bond screening effects. 
Ce// is highly heterogeneous in these molecules, and the 
largest values arise from the strongly overlapping acc 
and TTcc bonds. Consequently, t^ffipec) (1-9 and 2.6 ) 
and eeff{7rcc) ( 1 - 8 ) in C 2 H 2 and C 2 H 4 are significantly 
larger than eeff {cFcH) (1.4^ 1.5) and ee//(crcc) in C 2 i 76 
(1.5). 

For systems with a finite gap, it has been proved that 
(a) WRn(r) decays exponentially with |r —R| [23, 24]; (b) 


(a) (c) 


c-c 


C-H 




(r-R)/a 


FIG. 1: Locality of AWRR/(r) in an ethylene oligomer 

(Ci 9 iL 4 o). (a) Wannier orbitals of C-C and C-H bonds, (b) 
The exponential decay of AITRR/(r) as a function of |r — R'| 
(only shown for C-H bond) and (c) of 11 AWrr/ 11 as a function 
of |R — R'|. Solid lines in (c) indicate a linear fit. 



FIG. 2: The (local) density response of silicon under a uni¬ 
form electronic field in the [001] direction, (a) Contour plot 
of the local density response in the [110] plane of a Si-Si bond 
(inset: total density response), (b) Density response profile 
along the Si-Si bond. 


the Hamiltonian matrix SYin.Yi'n' decays exponentially 
with |R — R'l [25]; and (c) the density matrix (r|P^|r') 
decays exponentially with |r — r'| [25-27]. The local¬ 
ity of AWRn,R'n'(r) can be derived accordingly, (i) For 
a regular potential that is not exponentially divergent, 
PcAWH/R/n'(r) decays exponentially with |r — R'|. This 
can be easily understood from the locality of H7Rn(r) 
and Py, as Pc = / — Py (ii) [e — {H P aPy)I]^^, de¬ 
cays exponentially with |R —R'|, a direct consequence 
of the locality of ^Rn,R'n' [34]. (hi) Putting (i) and (ii) 
together, we conclude that ALFRn,R'n'(i*) decays expo¬ 
nentially with |r — R'l for a given R and R' pair, and 

its two-norm, | |AirR„,R /„/11 = ^ 

decays exponentially with |R — R'|. 

To validate these arguments, we considered an ethy¬ 
lene oligomer (C 19 P 40 with unit length a) containing two 
types of WFs, C-C and C-H bonds, as shown in Fig. la. 
To quantify the locality of AWRnR'n'? we consider a 
linear potential Vs{r) = x applied along the molecule. 
Both AWR„R/„4r) (Fig. lb) and ||AhFR„,R/„/|| (Fig. Ic) 
clearly exhibit an exponential decay. While the decay 
length of the former ( 0.6 ^ 1 . 0 a) is governed by the lo¬ 
cality of the WFs (0.3 ^ 0.5a) and the density matrix, 
the latter (0.5 ^ 0.7a) is controlled by the locality of 
%nR'n' (see Fig. S4 in SM). 

The strong spatial localization of the response density 
also exists in 3D crystals, as shown in bulk silicon in 
Fig. 2. Under a uniform electronic field, although the 
total response density is delocalized, 82% local density 
response along the [ill] direction is confined within one 
Si-Si bond (see Fig. 2 b). 
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An important outcome of the exponential localization 
of and AWun.R'n'ir) is that Xr„,r'„' is exponen¬ 

tially localized. It supports the physical picture of as a 
linear combination of coupled local response modes cen¬ 
tered on WFs. Consequently, efficient algorithms can be 
developed to compute taking advantage of its sparse 
representation in the real space. To prove this point, we 
construct a tight-binding model of for silicon with a lo¬ 
cal basis set, and compute its eigenvalue spectrum in the 
Brillouin zone using the dielectric band structure (DBS) 
interpolation, in analogy to the electronic band structure 
interpolation using Wannier functions [12]. To our best 
knowledge, this is the first demonstration of DBS inter¬ 
polation for a covalent bonded crystal. 

We first construct local basis set of denoted by 
as approximate eigenvectors of where 

hopping terms eRn,R'n' (Rri 7 ^ R'n') are switched off. 
This procedure decouples equations in Eq. 6 , making 
them easy to solve computationally. Then we calcu¬ 
late the hopping matrix X and overlap matrix S in 
real space, Vy'(R) = and 5 , 7 '(R) = 

Finally we Fourier transform them into mo- 
mentum space, 

X(q) = 

R 

5(q) = (10) 

R 

The interpolated DBS of x^ is the solution of the gener¬ 
alized eigenvalue problem in non-orthogonal basis: 

X{q)v = X{q)S{q)v. (11) 

We demonstrated the DBS interpolation method with 
bulk silicon using a 4 x 4 x 4 /c-mesh for the direct calcu¬ 
lation, and a 4 X 4 X 4 super cell with T-point sampling to 
construct the local basis set, X and S in real space [35]. 
Twenty five local eigenmodes per WF are sufficient to 
reach the numerical convergence of the first 25 bands of 
the total x^- The Wannier orbitals and the first ten di¬ 
electric basis functions of bulk silicon are shown in Fig. 3a 
and 3b. Excellent agreement (maximum absolute error: 
2.6 X was found between the direct calculation 

and interpolated DBS in Fig. 3c, which highlights the 
validity and accuracy of the tight-binding model of x^- 

The local basis set is the key to the efficient construc¬ 
tion of x^- Thanks to the exponential localization of 
XRn R'no Eq. 6 in principle can be solved within a sub¬ 
space containing Un neighboring Wannier orbitals of the 
local perturbation. Since is system size independent, 
the computational cost grows as 0{N^ln{N)) with N 
being the size of the system, which is a tremendous im¬ 
provement over the standard 0{N^) scaling (see SM). 
This quadratic scaling method is a promising starting 
point to develop low scaling algorithms for excited state 


(a) (b) 

^5® ®%^ ®^^ ®^4 ®^^ ®V« 

®^ ®^ 


®kA^ ®kA^ 

^ ^ it 

®^f^ ®^^ 



FIG. 3: DBS (x°) interpolation for bulk silicon using local 
basis, (a) The WF in bulk silicon, (b) The first 10 eigenmodes 
of the local response function, (c) Comparison of the DBS 
obtained from a direct calculation with interpolation using 
the lowest 25 basis functions per WF. 


problems, where the evaluation of x^ is the major bot¬ 
tleneck to construct, e.g., crpa = 1 — vx^ ot x [28-30]. 

The theoretical framework of local EDRFs can be ex¬ 
tended to metallic systems, as DFPT can treat metallic 
systems in general [18], but the decay rate of key quanti¬ 
ties can behave qualitatively differently from semiconduc¬ 
tors or insulators. In ID free electrons, the Wannier or¬ 
bital of the occupied portion of bands decays at r~^ [25]. 
For Wannier orbitals of disentangled bands, e.g., narrow 
transition metal d-bands, both numerical evidence and 
the analogy with the isolated composite case suggest the 
possibility of the exponential localization [13]. On the 
other hand, the decay of density matrix is expected to be 
algebraic at zero temperature, and exponential at finite 
temperature [26, 27]. The locality of EDRFs in metallic 
systems is therefore more subtle, and warrants further 
investigation. 

In conclusion, we present a local representation of 
EDRFs based on the concept of the dielectric response 
of WFs. This method allows us to perform fully ab ini¬ 
tio real space partition of EDRFs, and analyze excited 
state properties, e.g., polarizability, in terms of chem¬ 
ical bonds. In systems with a gap, we proved that the 
bare response function, Xr^ R'no decays exponentially in 
real space. This “near-sightness” is central to the physi¬ 
cal understanding of EDRFs and the development of low 
scaling algorithms for excited state problems. 

This work was performed at the Center for Functional 
Nanomaterials, which is a U.S. DOE Office of Science 
User Facility, at Brookhaven National Laboratory under 
Contract No. de-sc0012704. 
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ELECTRONIC STRUCTURE OF C2i7n, (n = 2,4,6) 

The polarizability decomposition of acetylene {C 2 H 2 ), ethylene {C 2 H 4 ), and ethane {C 2 HQ) was performed with 
WannierOO [1] and a modified version of Quantum ESPRESSO [2] using the Perdew-Zunger functional [3] in a 
35 X 18 X 18 Uq super-cell. The norm-conserving pseudo potential was used with energy cutoffs of 40 and 160 Ry 
for the wavefunction and charge density, respectively. The maximally localized Wannier functions (MLWFs) of these 
molecules are shown in Fig. SI. For C 2 H 2 and C 2 H 4 , a regular construction of WFs results in a mixture of tt and 
a bonds. Instead, a (from px) and tt (from py and Pz) orbitals were constructed separately, where x is along the 
direction of the C-C bond. 



•6 

(a) C 2 H 2 


(b) C2i74 (c) C 2 HQ 


FIG. SI: The WFs of C 2 H 2 (two ixcc and one acc bonds), C 2 HA (one ixcc and one acc bonds), and C 2 HQ (one acc bond). 


The magnitude of the unscreened bond polarizability (a^) is inversely related to the transition energy (AF) between 
the occupied states and the onset of the unoccupied manifold, i.e., the lowest unoccupied molecular orbital (LUMO). 
The relative energy of different chemical bonds is indicated in the projected density of states (PDOS) onto Wannier 
orbitals in Fig. S2. Since WFs are not eigenfunctions of the Kohn-Sham Hamiltonian, their PDOS can spread in an 
energy range, making it hard to compare between different WFs. Here we define the average energy of a WF as the 
on-site Hamiltonian matrix element, {Wn\H\Wn)^ indicated by the colored dashed lines in Fig. S2. Within the same 
molecule, in average AE{7rcc) < AE{CH) < AE{acc)^ which implies that a^{7rcc) > > a^(acc)‘ If also 

explains the increasing order of bond polarizabilities in C 2 H 2 , C 2 H 4 , and C 2 HQ for the same kind of bonds. 
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FIG. S2: Density of states projected onto WFs of C 2 H 2 , C 2 H 4 , and C 2 HQ. Colored dashed lines indicate the average energy 
levels of Wannier orbitals. 
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LOCALITY OF THE RESPONSE OF WANNIER FUNCTIONS 

The locality of AWRn,R'n'(r) has two meanings as described in the main text. First AWRn,R'n'(r) decays expo¬ 
nentially with |r — R'l for a given pair of R and R'. Second the two-norm, 11AlURn,R'n'11? decays exponentially with 
|R-R'|. 

In fact, the first locality is closely related to that of the WFs and the density matrix. As shown in Fig. S3, the 
quantity F(r) = PcAVsWRnir) exponentially decays in the real space, and the decay length (0.6a and 1.0a for C-H 
and C-C bonds, respectively) is approximately twice as its corresponding WF (0.3a and 0.5a for C-H and C-C bonds, 
respectively, not shown). This is because the projector, Pc(r, r') = I — P^(r, r'), has a similar decay length as the 
WF, therefore broadens the spread of the WF by two in P(r). Comparing Fig. S3 and Fig. lb in the main text, we 
observe that the locality of AH/RnR'n(i*) is similar to that of P(r), which means that the operator [s — {H — aPy)I]~^ 
has very narrow spread in the real space. 



(r-R)/a 

FIG. S3: The exponential decay of P(r), which results from the locality of both Wannier functions and the density matrix. 

On the other hand, the decay behavior of 11 AlFRn,R'n'11 dominated by the hopping term, £Rn,R'n'- We show the 
exponential decay of ejin,ii'n' in Fig. S4a, and the strong correlation between eun^ii'n' and 11 AlFRn,R'n'11 in Fig. S4b. 



FIG. S4: (a) The exponential decay of the Hamiltonian matrix in the Wannier basis, (b) The strong correlation between the 

magnitude of the local response and the hopping terms suggests that the rapid decay of the charge-flow is due to the decay of 
hopping terms between different WFs. 
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SCALING ANALYSIS OF 


In this section, we derive the numerical scaling to compute We will only consider the static limit, as the compu¬ 
tational cost to calculate imaginary frequency follows the same scaling. Here we focus on the computationally 

challenging case of non-periodic extended systems. The computational model normally contains a large number of 
electrons, A, in a super cell. For this reason, we consider only the F-point in the Brillouin zone, and drop k and 
q in the notation. Under the plane-wave implementation, suppose such a system has Ny occupied states and Ny 
unoccupied states at a given kinetic energy cutoff. The number of occupied Wannier orbital, Nyj, is the same as Ny. 
The size of the plane-wave grid is denoted by Nq- Here and in the following we use capital N for variables that scale 
with the system size, and use lower case n for variables independent of the system size. 

The most straightforward way to compute is to use the Adler-Wiser formula [4, 5]: 
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For semiconductors and insulators, Eq. 1 requires Ny x Ny Fast Fourier Transforms to calculate all pairs of co-density 
'0*(r)'0j(r) in the G space, whose cost is proportional to Ny x Ne X NgIuNg- Evaluating the full matrix of x^ requires 
four loops for G, G', i, and j, so the leading order of the scaling is 0{N^). Besides, the direct approach needs to store 
Nq elements. Both the computational and storage expenses become major bottlenecks for large systems. 

Alternatively, as illustrated in the main text, we could represent x^ i^i ^ local basis set, {|^f)}, made of approximate 
eigenmodes of local response. ) refers to the i-th approximate eigenvector of the local response of the n-th WF, 
If we need rim eigenmodes for each WF, the cost of the storage is proportional to rim x x Ng- Since rim is 
independent of the system size, the cost to store x^ grows as 0 {N^), but with a much smaller prefactor as compared 
to the direct evaluation of the Adler-Wiser formula. 

Now we analyze the computational cost of the local basis set approach. In general {|^f)} are obtained by iteratively 
solving the generalized Sterheimer equation and update AVg for several times (typical in the order of ten), until the 
new AVg is an eigenstate of x^- The matrix form of this equation reads: 


/ ^ 0,0 — H ... 




( AWo \ / PcAU.Wo 


\ ••• - H J \ AW 


V Pc^VgWn 


( 2 ) 


We construct the local basis set by removing the off-diagonal elements of e, and solve the eigenvectors of 


{ei,i-H)AWi = PcAVgWi. (3) 

To apply projector = I — | W)(W| in principle requires Ny x Ng operations, however, we can benefit from the 
exponential locality of W? and truncate the summation of the Wannier centers to neighbors within a properly 
chosen cutoff distance of the local perturbation without losing accuracy. The computational cost to calculate the 
right hand side of Eq. 3 therefore scales as x Ng, and the iterative solution of AW scales as Ng In The cost 
to construct the full local basis set will acquire a prefactor of rim x which gives an overall 0 {N^) scaling. 

To compute Xf'P = |x^|^7 ^ nsed in the dielectric band structure interpolation, one has to solve Nyj coupled 

equations in Eq. 2. Again, thanks to the exponential localization of x^, only rin coupled equations needs to be solved 
for each |^^ )• Therefore the leading cost to construct X is times of the cost of building the basis set, which again 
grows as O(A^). Compared to the O(A^) scaling of the standard Adler-Wiser formula, the local basis set approach 
has a tremendous saving in both computer time and memory. 
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